\documentclass[%
 prd,
%jmp,
 floatfix,
 amsmath,amssymb,
 preprint,%
 %reprint,%
%author-year,%
%author-numerical,%
]{revtex4-1}

\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math

\begin{document}

\title{Moments and Polarisabilities}

\author{Thomas Primer}
\affiliation{CSSM, University of Adelaide, SA 5005, Australia}

\author{Waseem Kamleh}
\affiliation{CSSM, University of Adelaide, SA 5005, Australia}

\author{Derek Leinweber}
\affiliation{CSSM, University of Adelaide, SA 5005, Australia}

\date{\today}% It is always \today, today,
             %  but any date may be explicitly specified

\begin{abstract}
We present results on the magnetic moment and magnetic polarisability of the neutron.
These results are calculated using the background field method on $32^3\times 64$ dynamical QCD lattices provided by the PACS-CS collaboration as part of the ILDG.
We use a uniform background magnetic field quantised by the periodic spatial volume.
We investigate ways to improve the effective mass plots used to calculate magnetic polarisabilities, including the use of correlation matrix techniques with various source smearings.
\end{abstract}

\pacs{Valid PACS appear here}% PACS, the Physics and Astronomy
                             % Classification Scheme.
\keywords{Suggested keywords}%Use showkeys class option if keyword
                              %display desired

\maketitle

\section{\label{intro}Introduction}

The magnetic moment and magnetic polarisability are fundamental properties of a particle that describe its response to a static external magnetic field.
The ability to calculate these properties is a useful test of lattice QCD.
There are two well known techniques for calculating magnetic moments on the lattice.
One is to use three-point functions to get electromagnetic form factors which can then be converted into magnetic moments \citep{Boinepalli:2006}\citep{QCDSF:2011}.
The other is the background field method, which uses a phase factor on the gauge links to induce a constant external field across the whole lattice \citep{Bernard:1982}\citep{Smit:1987}\citep{Burkardt:1996}\citep{Zhou:2003}\citep{Lee:2005}\citep{Lee:2006}\citep{Lee:2008}.
This external field causes a mass shift from which the magnetic moment and polarisability can be derived by making use of  the following energy-field relation,
	\begin{equation}\label{erelation}
	E(B) = M_N + \vec\mu \cdot \vec{B} + \frac{e\vert{B}\vert}{2M_N} - \frac{4\pi}{2}\beta B^2 + \mathcal{O}(B^3).
	\end{equation}

%When deriving the background field method, the need to have a consistent magnetic field across the periodic boundary leads to a quantisation condition which limits the available choices of magnetic field strength based on the size of the lattice.
When deriving the background field method on a periodic lattice there arises a quantisation condition which limits the available choices of magnetic field strength based on the size of the lattice.
If the lattice is too small the field will be large and higher order terms in the energy relation (Eq.~\eqref{erelation}) will begin to dominate.
Previous calculations have avoided this problem by using Dirichlet boundary conditions and a linearised form of the phase factor, which allows for an arbitrary choice of field strength.
Others have used the exponential phase, but instead of correcting the value of the field at the boundary they put the quark origin at the centre of the lattice and hope that the boundary is far enough away for the effects of the discontinuity to be small.
Using either of these methods introduces finite volume errors which can be hard to predict.
Our calculation is the first to use periodic boundary conditions and the quantised exponential phase factor.
We present results for both the magnetic moment and magnetic polarisability of the neutron.


\section{\label{method}Background Field Method}

We make use of the background field method to simulate a constant magnetic field along one axis.
The technique is formulated on the lattice by first considering the continuum case, where the covariant derivative is modified by the addition of a minimal electromagnetic coupling
\begin{equation}
D_{\mu} = \partial_{\mu}+gG_{\mu}+qA_{\mu},
\end{equation}
where $A_{\mu}$ is the electromagnetic four-potential and $q$ is the charge on the fermion field.
On the lattice this is equivalent to multiplying the usual gauge links by a simple phase factor
\begin{equation}
U_{\mu}^{(B)}(x)=\exp(iaqA_{\mu}(x)).
\end{equation}

To obtain a uniform magnetic field along the z-axis we note that $\vec{B} = \vec{\nabla} \times \vec{A},$ and hence 
\begin{equation}\label{Bz}
B_z = \partial_x A_y - \partial_y A_x.
\end{equation}
Note that this equation does not specify the gauge potential uniquely.
There are multiple valid choices of $A_\mu$ that give rise to the same field, however they are equivalent up to a gauge transformation.
We choose $A_x = By$ to produce a constant magnetic field of magnitude $B$ in the $z$ direction.
Examining a single plaquette in the $(\mu,\nu) = (x,y)$ plane shows that this gives the desired field for a general point on the lattice.
However on a finite lattice $(0 \le x \le N_x-1), (0 \le y \le N_y-1)$ there is a discontinuity at the boundary due to the periodic boundary
conditions.
In order to fix this problem we make use of the $\partial_x A_y$ term from equation \eqref{Bz}, giving $A_y$ the following values,
\begin{equation}
A_y(x,y) = \begin{cases} 0 & \mathrm{for\ } y < N_y-1 \\
	           -N_y Bx & \mathrm{for\ } y = N_y-1, \end{cases}.
\end{equation}
This ensures that we now get the required value at the $y=N_y-1$ boundary.

There is then the issue of the double boundary, $x = N_x-1$ and $y=N_y-1$, where the plaquette only has the required value under the condition
$\exp(-ia^2 qBN_x N_y)=1$.
This gives rise to the quantisation condition which limits the choices of magnetic field strength based on the lattice size
\begin{equation} \label{quantisationcondition}
a^2 qB = \frac{2\pi n}{N_x N_y},
\end{equation}
where $n$ is an integer specifying the field strength in multiples of the minimum field strength quantum.


\section{\label{details}Simulation Details}

These calculations have been performed on 2+1 flavour dynamical-fermion configurations provided by the PACS-CS group  \citep{PACS-CS} through the ILDG \citep{ILDG}.
These are $32^3 \times 64$ lattices using a clover fermion action and Iwasaki gauge action with $\beta=1.9$ and physical lattice spacing $a=0.0907(13)$ fm.
We use four values of the light quark hopping parameter, $\kappa_{ud} =$ 0.13700, 0.13727, 0.13754, 0.13770, corresponding to the pion masses $m_\pi $= 702, 572, 413, 293 MeV.
[We also have gauge configurations for $\kappa_{ud} =$ 0.13781, however calculating background field propagators at this mass has proved difficult. The light mass already leads to the matrix inversions taking a long time and the addition of the background field means that many of them have difficulty converging at all.]

In order to get correlation functions at four different magnetic field strengths we calculated propagators at six non-zero field strengths, $qBa^2 = $+0.0061, -0.0123, +0.0184, +0.0245, -0.0368, -0.0492.
These were then combined using a down quark propagator with field $-B$ and up quark propagator with field $2B$ to get an overall field of $3B$.
Unless specified otherwise we used the interpolating field $\chi_1 = (u^TC\gamma_5d)u$ with 100 sweeps of stout smearing at the source.

It should be noted that the configurations are dynamical only in the QCD sense, there was no magnetic field included when they were calculated.
The background field can be put on the sea quarks by performing a separate calculation for each field strength, but this is obviously very computationally expensive.
It also destroys the correlations between the different field strengths which would lead to much larger errors in the energy shifts used to calculate moments and polarisabilities.
There also exist techniques for a reweighting of configurations in order to correct for the background field \cite{Freeman}, however these have not been employed in this work.
It is believed that the corrections due to the effect of the background field on the sea quarks will be small.

\section{\label{moment}Magnetic Moment}

When a charge or system of charges with angular momentum is placed in an external magnetic field it is energetically favourable to have its axis of rotation either aligned or anti-aligned with the direction of the field.
The tendency of the system to align with the field is proportional to the magnetic moment of the system and the strength of the field.
The ground state of a hadron like the neutron has no orbital angular momentum, so the magnetic moment is due only to the intrinsic spin angular momentum of its valence quarks and sea quarks.

We calculate correlation functions containing spin-up and spin-down components in the (1,1) and (2,2) positions of the dirac matrix respectively.
For a magnetic field aligned to the axis of the spin we see the magnetic moment manifested as a shift in the energy which has the same magnitude, but opposite sign, for spin-up and spin-down.

%\subsection{Method}
We make use of the sign difference in the energy shift between spin-up and spin-down in order to isolate the magnetic moment term from the expansion of the energy, taking the difference of the spins,
\[ \frac{1}{2}\left(E_\uparrow(B) - E_\downarrow(B)\right) = \mu B. \]
In terms of correlation functions there are multiple valid ways of achieving this, for example fitting the energy and then taking the difference or combining correlation functions and then fitting.
By taking a combination of correlation functions before fitting for the energy the statistical error is greatly reduced.
This is because the errors are highly correlated between the zero and non-zero field correlation functions, meaning the fluctuations do not change significantly due to the field.
The combination required for isolating the moment term can be written as,
\begin{equation}
\delta E(B) = \frac{1}{2}\left(\ln\left(\frac{G_{\uparrow}(B,t)}{G_{\uparrow}(0,t)}\frac{G_{\downarrow}(0,t)}{G_{\downarrow}(B,t)}\right)\right)_{\rm fit}.
\end{equation}
The inclusion of the bare correlation functions without a magnetic field in this expression is not strictly necessary, but it is useful in correcting for the small statistical difference between spin-up and spin-down zero field energies and making the zero field point zero by construction.
We also redefine spin-up to mean aligned with the magnetic field and spin-down to mean anti-aligned to the field so that we can treat all the fields as positive.


\subsection{Results}

\begin{figure}%[t]\hspace{-24pt}
\includegraphics[width=0.35\textwidth,angle=90]{Kappa1Diff.pdf}
\caption{(Colour online) Spin-differenced effective mass for the heaviest quark mass at all four field strengths.}
\label{diffplots}
\end{figure}

Figure \ref{diffplots} shows the effective mass shift from the difference of spin-up and spin-down for the heaviest quark mass at all four non-zero magnetic field strengths.
These generally display excellent plateaus, with similar results for the other quark masses.
At the two higher field strengths there is a small drift from the early plateau at 18-20 to a second plateau at 23-24.
This leads to a slight difference in the value of the mass shift depending on the choice of fit window, however this has little effect on the magnetic moment result as we will now demonstrate.

\begin{figure}%[t]\hspace{-46pt}
\includegraphics[width=0.35\textwidth,angle=90]{NeutronDiffPlot.pdf}
\caption{(Colour online) Fits of the spin-difference mass shift to the field strength at each quark mass. The solid line is a purely linear fit to just the first two points and the dashed line is a linear plus cubic fit to all four points.}
\label{difffits}
\end{figure}

Figure \ref{difffits} shows the spin-difference effective mass shifts plotted against the magnetic field strength.
These are fit to a linear coefficient which gives the magnetic moment.
In order to fit the largest field strength, and to a lesser extent the second largest, we had to include a cubic term in the fit.
The cubic term is able to absorb some variation in the mass shifts at the higher field strengths, which is why the shifting plateaus in Fig.~\ref{diffplots} don't really matter.
This means that the first two points are the main drivers of the linear coefficient and therefore the magnetic moment.
We also performed a purely linear fit to only the first two points and found that the linear coefficients agreed well within errors.
Since the fit is naturally constrained to go through zero, two non-zero field points are enough to give us confidence that our field strengths are small enough to make the higher order contributions negligible.
Only at the lightest mass does this begin to maybe come into question, potentially leading to a slight underestimate of the magnetic moment value.

\begin{figure}%[t]\hspace{-24pt}
\includegraphics[width=0.35\textwidth,angle=90]{MomentWithExtrapolation.pdf}
\caption{(Colour online) Neutron magnetic moment as a function of pion mass squared. The left most point gives the experimental value \cite{PDG:2012}. The dashed line is a fit of the dynamical points to $\mu_n = \mu_0/(1-\frac{\chi_n}{\mu_0} m_\pi + c m_\pi^2)$.} %$\mu=a+bm_\pi$.}
\label{momentplot}
\end{figure}

Figure \ref{momentplot} shows the neutron magnetic moment results, compared with a previous quenched calculation. 
The two sets of results agree within errors.
The line is a chiral fit of the form,
\[ \mu_n(m_\pi) = \mu_0/(1-\frac{\chi_n}{\mu_0} m_\pi + c m_\pi^2) \]
%\[ \mu=a+bm_\pi, \]
where $\chi_n = \beta_n^{(\pi)}(m_N/8\pi f^2)$ has the value $4.41$, calculated from one-loop corrected estimates of $\beta_n^{(\pi)}$ in \citep{HackettJones:2000}\citep{Leinweber:1999}, and $\mu_0$ and $c$ are fit freely.
This results in an extrapolated value at the physical limit of $-1.82\, \mu_N$, only slightly above the physical value of $-1.91\, \mu_N$.
We also expect there to be significant finite volume effects due to the small $m_\pi L$, but those have not been examined here.
The magnetic moment is given in units of nuclear magnetons ($\mu_N = \frac{e\hbar}{2M_N}$) which is achieved by scaling the result by the physical nucleon mass.% [check for more precise way to say this?]. [needs some sort of closing sentence for this section I think].

\section{\label{polarisability}Magnetic Polarisability}

The magnetic polarisability is a measure of the deformation of a non-point particle when it is placed in a magnetic field.
This deformation causes a change in the energy which we can measure using the background field method.
The effect of the magnetic polarisability is second order in $B$.
This means that at the "small" field strengths we are using the effect is much smaller than that due to the magnetic moment, which can make it hard to measure.
It also makes it more important to use the full exponential phase factor, since the errors introduced by the linearised form are also at order $B^2$\cite{Guerrero:2008}.

%\subsection{Method}

To extract the polarisability from the energy we take the average of spin-up and spin-down to remove the magnetic moment term and explicity subtract the zero-field mass.
\[ \frac{1}{2}\big(E_\uparrow(B)-E_\uparrow(0)) + (E_\downarrow(B)-E_\downarrow(0))\big) = \frac{e\vert{B}\vert}{2M_N}- \frac{4\pi}{2}\beta B^2. \]
This leaves us with the polarisability term, but also with another term due to the Landau energy.
This energy arises from the quantisation of orbits for charged particles in magnetic fields and can't be isolated from the relevant polarisability term.
This term makes it difficult to calculate magnetic polarisabilities of charged particles because there is not only the ground state Landau energy but also a tower of Landau levels.
The need for small field strengths makes the Landau level problem even worse because it means the Landau levels are closer together, which makes it take longer in euclidean time for the levels above the ground state to be exponentially suppressed \citep{Tiburzi:2012}.
Fortunately for a neutral particle like the neutron the Landau term is zero and can be ignored.

As with the magnetic moment we construct ratios of correlation functions which we then fit for an effective mass,
\begin{equation}\label{polratio}
\delta E(B) = \frac{1}{2}\left(\ln\left(\frac{G_{\uparrow}(B,t)}{G_{\uparrow}(0,t)}\frac{G_{\downarrow}(B,t)}{G_{\downarrow}(0,t)}\right)\right)_{\rm fit}.
\end{equation}
In this case the zero-field correlators are necessary to remove the bare neutron mass.
Combining the correlation functions before fitting is especially important in the polarisbility case because the mass shift is smaller than the errors on the zero-field mass.
This means that if the correlated errors were not allowed to cancel before the fit we would not get any signal for the shift due to the polarisability. [show plot of this maybe?]\\

\subsection{Results}

\begin{figure}%[t]\hspace{-42pt}
\includegraphics[width=0.35\textwidth,angle=90]{Kappa1Sum.pdf}
\caption{(Colour online) Spin-averaged effective mass shift for the heaviest quark mass at all four field strengths.}
\label{sumplots}
\end{figure}

Figure \ref{sumplots} gives the spin-averaged effective mass shift for the heaviest quark mass.
Unlike the spin-difference case the plateau behaviour is quite poor, with a fairly constant downward slope that begins to plateau only when significant noise is appearing.
The situation is similarly poor at other quark masses.
The plot shows that at each field strength the mass shift starts at approximately zero and grows with Euclidean time.
This is because the polarised state has a deformed wavefunction which has poor overlap with our spherically symmetric smeared source.
It takes time for the presence of the background field to deform the particle to its fully polarised state [use plot of Dale's work showing spherical wavefunction at time zero which then squishes?].

Since the late plateaus are related to the shape of the wavefunction we believed we could potentially improve our results by finding a better suited source.
We experimented with a number of different source smearings, trying 16 and 35 sweeps in addition to our usual 100.
We also tried a point source on the basis that it should have no bias towards any shape and may therefore reach the required form more quickly.
Figure \ref{upcompare} shows the mass shift due to the field at all smearings for the heaviest quark mass at the smallest field.
The three smeared sources have different asymptotic behaviour but agree well within errors by time slice 25, just before the signal is lost to noise.
The point source has large excited state contributions and does not match the other sources before losing signal.
This suggests that although the excited states cancel in the spin-average to give decent looking plateaus they cannot be relied on because it is not reaching the ground state.

We notice from these plots that the best plateau behaviour comes from 16 sweeps of smearing for spin-up and 100 sweeps of smearing for spin-down.
To take advantage of this we constructed a spin-average from the spin-up correlation funciton with 16 sweeps and the spin-down correlation function with 100 sweeps.
We found that this did in fact result in improved plateau behaviour, as shown in Figure \ref{sumcompare}, leading us to investigate further the possibility of combining different source smearings.

The variational method \cite{Mahbub:2011} involves using an $n \times n$ correlation matrix $G_{ij}(t)$ constructed from different sources and sinks to solve a pair of eigenvalue equations.
The eigenvectors $u_j^\alpha$ and $v_i^\alpha$ can then by used to project out energy eigenstates $\alpha$ to effectively isolate the $n-1$ lowest energy states.
\begin{equation}\label{eigenvector}
G^\alpha(t) = v_i^\alpha G_{ij}(t)u_j^\alpha
\end{equation}
Combining correlation matrix techniques with the background field method introduces new considerations.
Normally the eigenvector analysis is performed on spin-averaged correlation functions because spin-up and spin-down are the same up to statistics.
For a background field calculation we need to consider spin-up and spin-down separately, with each field strength getting its own eigenvector equation to solve.
We then construct the same ratio as in Eq.~\ref{polratio} but using the projected correltion functions from Eq.~\ref{eigenvector}.

Performing a variational analysis using a $2\times 2$ correlation matrix made from 16 and 100 sweeps of smearing allows us to get the best possible combination of each for isolating the ground state.
When we did so we found that the plateau behaviour was not noticably better than using either of the smearings alone.
This suggests that the good plateau in Figure \ref{sumcompare} was simply a conspiracy of conveniently cancelling excited states that only looked good by coincidence.
We also tried other various combinations of sources with different smearings and interpolating fields in the correlation matrix and found that none gave a statistically significant improvement.

\begin{figure}[t!]%[t]\hspace{-46pt}
\includegraphics[width=0.35\textwidth,angle=90]{SmearUpCompare2.pdf}\vspace{-12pt}\\%\hspace{-46pt}
\includegraphics[width=0.35\textwidth,angle=90]{SmearDownCompare2.pdf}%\hspace{-24pt}
\caption{(Colour online) Mass shift at the smallest field strength. For spin-up we have from top to bottom we have 100, 35, 16 sweeps of smearing and point source, for spin-down the order is reversed.}
\label{upcompare}
\end{figure}%\hspace{-46pt}

\begin{figure}%[t]\hspace{-42pt}
\includegraphics[width=0.35\textwidth,angle=90]{SumCompareCombine.pdf}%\hspace{24pt}
\caption{(Colour online) Spin-averaged effective mass shift for the heaviest quark mass at the smallest background field using four different source; Point source, 16, 35 and 100 sweeps of smearing.}
\label{sumcompare}
\end{figure}

%\begin{figure}[t]\hspace{-46pt}
%\includegraphics[width=0.4\textwidth,angle=90]{SmearDownCompare2.pdf}%\hspace{24pt}
%\caption{(Colour online) Spin-down mass shift at the smallest field strength. From top to bottom we have point source, 16, 35 and 100 sweeps of smearing.}
%\label{downcompare}
%\end{figure}

Since we could not achieve good plateaus for fitting, we extracted energy shifts by simply taking the value at the point just before noise starts to dominate the signal, where we are confident the excited states have been suppressed because the all different smearings agree.
Figure \ref{sumfits} shows fits of these spin-averaged mass shifts to the field strength.
In addition to the quadratic polarisability term we required a quartic term in order to fit the higher field strengths.
This higher order term is fairly negligible at the heaviest three masses but starts to become significant at the fourth mass.

\begin{figure}%[t]\hspace{-46pt}
\includegraphics[width=0.35\textwidth,angle=90]{SumPlotsNew.pdf}%\hspace{24pt}
\caption{(Colour online) Fits of the spin-averaged mass shifts to the background field strength. The dashed line is a fit of all the point to a quadratic plus quartic.}
\label{sumfits}
\end{figure}

Figure \ref{polarisabilityplot} shows our neutron magnetic polarisability results compared with previous quenched results.
The quenched and dynamical results agree well within errors.
The dashed line shows a fit of the dynamical results to
\[ \beta = a + b m_\pi + c/m_\pi. \] % \frac{c}{m_\pi} \]
The extrapolated value is right in the middle of the experimental result, which has large errors bars.
With an improved calculation at near physical pion mass we could soon expect a lattice result that improves on the experimental value for the neutron magnetic polarisability.

In order to get better polarisability results in the future we believe the best approach would be to use anisotropic sources with their smearing tuned to the shape of the polarised particle.
This would allow the spin-averaged mass shift to start at near its final value and reach a plateau before statistical noise becomes too large.
The other approach would be to just massively increase the statistics in order to stave off noise long enough for clear plateaus to appear, however this is beyond our resources.

\begin{figure}[h]%[b]\hspace{-46pt}
\vspace{18pt}
\includegraphics[width=0.35\textwidth,angle=90]{PolarisabilityTest.pdf}%\hspace{24pt}
\caption{(Colour online) Neutron magnetic polarisability vs pion mass. The red point gives the experimental value \cite{PDG:2012}.The line represents a fit of the dynamical points to $\beta = a + bm_\pi + c/m_\pi$.}
\label{polarisabilityplot}
\end{figure}



\section{Conculsion}

We have performed the first calculations of the magnetic moment and magnetic polarsability of the neutron in a uniform background field on dynamical lattices.
Results for the magnetic moment are very clear and agree well with previous calculations.
The chiral extrapolation gives a value of $-1.82\, \mu_N$ at the physical pion mass.
Finite volume corrections will be required to get a final lattice value and will be made easier by results at lighter masses.

Magnetic polarisability calculations have proved more difficult due to late appearing plateaus.
We have calculated results which agree with previous quenched results and have an excellent approach to the experimental value.
We believe these polarisability calculations can be improved through the use of anisotropicly smeared sources.

\bibliography{Bibliography}{}
\bibliographystyle{ieeetr}

\end{document}
